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Abstract 

Using molecular dynamics simulations, we study the effect of pressure on the binding propensity 
of small RNAs by calculating the free energy barrier corresponding to the looped out conformations 
of unmatched base, which presumably acts as the binding sites for ligands. We find that the 
free energy associated with base looping out increases monotonically as the pressure is increased. 
Furthermore, we calculate the mean first passage time of conformational looping out of the base 
bulge using the diffusion of reaction coordinate associated with the base flipping on the underlying 
free energy surface. We find that the mean first passage time associated with bulge looping out 
increases slowly upon increasing pressures P upto 2 kbar but changes dramatically for P > 2 kbar. 
Finally, we discuss our results in the light of the role of hydration shell of water around RNA. 
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I. INTRODUCTION 



RNA molecules are very diverse both structurally and functionally [1]. Apart from having 
the regular helical purine-pyramidene base pairs, RNA molecules are also found to have many 
other secondary structures (motifs) such as loops, knots and bulges etc [2|, yj. The presence 
of such structural motifs is found to play a role in binding of different molecules to RNA (4]. 
For many protein binding RNAs, it was found that the frequency of adenosine bulge at the 
binding site is very high. The presence of a bulge may change the conformational flexibility 
of an RNA 5|, |6[ and hence more internal surface area of RNA is available for any chemistry. 
Moreover, the presence of the bulges does not only change the conformational flexibility but 
the bulges themselves may just flip out exposing the internal regions of a RNA to solvent 
and ligands. It has been shown that the bulge base looping out is highly sensitive to the 



bulge bases and their neighbors {5, 6], which makes the question of generality of any picture 
of base bulge looping out difficult. 

Recent state of the art computer simulations of small RNAs have shed light on the base 
bulge looping out processes {7- 9] in few of the possible cases. Specifically, these works have 
studied the free energy barriers associated with bulge base looping out process |7j ■ For exam- 
ple A. Barthel and M. Zacharias studied single uridine and adenosine bulge structures and 
their looping out using the torsional angle which measures the degree of looping out of bulge 
bases from the local helical plane. They find that the conformational free energy change in 
the case of adenosine bulge in a complete looping out process is higher by 1.5 kcal.mol -1 as 
compared to the adenosine bulge [7f, suggesting that in a base nonspecific binding process 
a structure with single uracil bulge base would have higher propensity to flip out of helical 
plane. Although, a wealth of literature is available on the base looping from the helical plane 
at ambient conditions; the changes in the kinetics of bulge base flipping is rather unexplored 
at conditions away from ambient conditions. 

It is known that water's hydrogen bond network and so the local structure of liquid water 
changes upon changing thermodynamic conditions, giving rise to anomalous changes in the 
dynamics and thermodynamics of water and aqueous systems. The structural stability and 
kinetics of proteins (where structural stability usually implies kinetically functional) as a 
function of pressure and temperature. The solvation barrier in the case of proteins plays an 
important role both in the hydrophobic collapse of the polypeptides as well as the stability of 



these structures as a function of pressure and temperature. In this paper, we will investigate 
another class of molecules which are involved in ligand binding kientics and are very stable 
at high pressures and temperatures. One can ask questions whether the structural stability 
in these cases will also imply kinetic functionality. We specifically study the binding of a 



purine A based molecules to small RNAs with a single A-bulge structural moti: 



. This study 



is relevant to understand the self-aminoacylation of small ribozymes [lOj, [ll| in different 
thermodynamic conditions. It is widely believed that for the purine A-based molecules such 
as ATP to bind to such a bulge, the A-bulge base has to flip out, which allows the purine 
A-based molecules to come and stack into the bulge configuration. In other words, for the 
binding to proceed the A-bulge has to overcome both the bending regidity and solvation 



energy barrier to so 
studied problem 
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vate in water. The solvation of different substances in water is a widely 
- Il7j . Indeed, studies of apolar solutes in water shows an ellipitic region 
in the pressure-temperature plane in which water behaves as a bad solvent and hence less 
solubility of these subtances [181]. However, study of solvation barrier of base flipping in case 
of RNA/DNA is not studied. 

In this paper, we study the effect of pressure on the kinetics of bulge base looping out 
of a double strand RNA with single adenosine bulge. In section II, we discuss the method, 
in section III we present the results for the free energy barrier associated with torsional 
deviation of the bulge base from the local helical backbone, in section IV we present a mean 
first passage time calculation associated with looping out process and finally we discuss and 
summarize the results. 



II. SYSTEM AND METHOD 



The energy minimized starting structure of the RNA with single A (5'-GGGGAGG-3'/5'- 
CCCCCC-3') bulge was created using the NAB/Nucgen module of the AmberlO program 
suite (Fig. H]). The RNA structure was then solvated in 4000 TIP3P water molecules such 
that there was about 1.0 nm space left between the boundary of the box and the RNA solute 
atoms. For electro-neutrality 14 Na+ counter ions were added to the solution. Energy min- 



imizations were carried using the steepest descent (1000 steps) in GROMACS3.3.3 19l. |20| 



with keeping the RNA atoms fixed. Simulations were carried out using the GROMACS3.3.3 
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FIG. 1: (a) A 2-dimensional schematics of the double strand RNA structure used in this study with 
a single A bulge in 5'-GGGGAGG-3'5'-CCCCCC-3'. (b) The three-dimensional structure and the 
Definition of the torsional angle (Cl'Cl'Cl'Nl) chosen as the reaction coordinate for calculation 
of the free energy. 

program with the Amber99 force field with a periodic boundary conditions and integra- 
tion time step of 2 fs. The long range electrostatic interactions were treated with the 
particle-mesh-ewald (PME) method. After minimization the system was slowly heated to 
temperature T=300 K with positional restraints. 

After the position restrained simulations, unrestrained Molecular Dynamics (MD) was 
carried out at four different values of the pressures P = 1,1000,2000,3000 atm. Thermal 
equilibrium at a constant temperature T = 300 K and different pressures was achieved using 
Berendsen thermostat and barostat respectively. The final equilibrated conformation was 
then used for as the starting conformation for umbrella sampling at different pressures. 
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To quantify the relative propensity of binding of RNA at different pressures, we chose 
the dihedral angle 9 (Cl'Cl'Cl'Nl) (see Fig. [1]) as the reaction coordinate for calculation of 
free energy. Since the conformational changes are very slow and an equilibilirium sampling 
of torsional angles require much larger time scales than computationally feasible, we use 
biasing potential to calculate the free energy. Umbrella sampling method was used to cal- 
culate the relative free energy of bulge base looping out conformations for the RNA shown 
in Fig. [TJ Harmonic umbrella biasing potentials Uuas,e l = k{9 — 9i) 2 with a force constant 
k = 0.05kcal.mol~ l .deg~ 2 were distributed uniformly along the reaction coordinate 9 at 
an interval A9 re f = 5°. Consecutive sampling windows of 9 were started from equilibrium 
structure of last run. For each values of the umbrella sampling window, we run a 2 ns sim- 
ulation and record the value of 9 every 0.2 ps. The final potential of mean force (PMF) was 



calculated using the WHAM (weighted histogram method) 21]. The unbiased probability 



distribution P(9) at a given temperature T under WHAM is given by 

p(0) = Ei=i w i(0) / j\ 

Y^-T n i e[F *~ Ubias,im/kBT 

where N sim is the number of sampling window (simulations), rii is the number of counts in 
the bin associated with 9, Uu a s,i is the biasing potential, and Fj free energy from simulation 
% and is given by 

Fi = -k B Tln[J2 P(9)e i ' Ub — i9)/kBT) } (2) 

®bins 

where 9u ns is the number of bins for the individual sampling window. The equations Q] and 
[2] are iterated to obtain the self consistent value of P{9). The value of P(9) depends on the 
time scale of simulations and hence long simulations are needed for a good convergence of 
the free energy. 

In Fig. [21 we show typical histograms along the reaction coordinate 9 for biasing potential 
centered at different values of 9 with A9 = 5°. We obtain smooth histograms for all the 
sampling windows for all the trajectories. 
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FIG. 2: A typical histogram for biasing potential centered at different values of 9. 

III. KINETICS OF SINGLE A BULGE AT AMBIENT AND ELEVATED PRES- 
SURES 

In Fig. |3l we show F{9) as a function of 9 for the RNA structure with A-bulge at 
P = 1 atm. Negative value of 9 corresponds to deviation towards minor groove while 
the positive values corresponds to deviation towards major groove. We find that F{9) has 
characteristic two minima centered around 9 ~ —20° and 9 ~ 30° as reported in earlier 
studies of single A bulge Note that the definition of the torsional angle 9 is different 
from the one used Ref. and hence different values of 9. As we can see from figure [3j the 
orientation of the A-bulge at more stable minimum is tilted slightly along the major groove 
while the second minimum at 9 ~ —20° is presumably due to the base triplet formation with 
the neighboring bases. The free energy difference between these two minima is ~ 3kcal/mol. 
suggesting that although 9 = 30° is relatively a more stable minimum configuration, the 
thermal fluctuations at T = 300 K 0.60 kcal/mol) is sufficient enough for the bulge to 
get displaced of the free energy minimum configuration. Due to a large free energy barrier 
(AF ~ 5kcal/mol) a complete looped out conformation of single A-bulge is less favorable 
and hence consistent with experimental fact that a single A does not bind efficiently with 
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the incoming ligands. 
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FIG. 3: Free energy as a function of the torsional angle 9 for single A bulge at P = 1 atm. The 
free energy profile associated with bulge flipping has two minima as found in earlier studies Q], 
separated by an energy barrier of 3.0 kcal/mol. 

Figure H] shows F{6) as a function of 9 for different P. We find that for pressures up 
to 2000 atm F(8) has the characteristic two minima as we find in the case of P = 1 atm. 
However, as the pressure in increased, the minimum at 9 m —20° becomes shallower and 
disappears for P > 2000 atm, suggesting that at P > 2000 atm the base triplet formation 
of the bulge base with the neighboring bases does not occur during the looping process. 
Moreover, the free energy barrier between the two minima changes just a little upon in- 
creasing pressure for P < 2000 atm. For P = 3000 atm, free energy barrier for the flipped 
out state changes drastically where AF m 8.7 kcal/mol, suggesting that the propensity of 
single A-bulge base looping out from the local helical plane would decrease upon increasing 
pressure and so the binding propensity of of incoming ligands. Furthermore, we note that 
although the free energy barrier associated with complete looped out conformations is finite 
is for all the pressures the final looped out state is not a free energy minimum state hence 
a RNA structure with single A bulge would not favor ligand binding. To verify this we ran 
20 ns long simulations to see if there is any bulge looping out event at different pressures 
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and we do not observe any. 




FIG. 4: Free energy as a function of the torsional angle 9 for different pressures. The free energy 
barrier between the two minima slowly disappears at pressure P > 2000; suggesting that the 
transient barrier that is produced at atmospheric pressures due to the partial base triplet hydrogen 
bonding of the bulge base with the neighboring bases is almost broken at higher pressures. 

In order to calculate the effective rate k P f f of transition from a stacked to looped out 



conformation, we use Langevin equation 
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23|. Assuming the diffusion of the reaction 



coordinate 9 on an underlying free energy surface the dynamics of 9 is governed by 

= -T^d e F(9) + f(t) (3) 
k B T 

where 9 is the torsional deviation and D is the diffusion constant and f(t) is the thermal 
noise with zero mean, (f(t)) = and delta function correlation, (f(t)f(0)) = 2D5(t). In 
the high friction limit, the probability p(9,t) of finding the system with reaction coordinate 
9 after time t is given by the Smoluchowsky equation: 

p(9,t) = Lp(9,t) (4) 

where L is the Fokker-Planck operator given by L = doe~ l3F ^ Ddoe l3F ^ and (3 = 1/ksT. 
The mean first passage time t{9i) associated with crossing the barrier from any coordinate 
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FIG. 5: Schematic of free energy as a function of the torsional angle 9 used for the calculation of 
the effective time scale of bulge base flipping kinetics. 

6i to final state 9f is given by (see Fig. [5]) 

t(9i) = f f dy±- f dx eOW)-^))), (5) 

JOi D Jg r 

where 9 r and 9 S denote the reflecting and absorbing boundaries respectively. We choose 
9f — 80° as the final looped out state and initial state is chosen to the values of 9 where the 
free energy curve has the deepest minimum for respective pressures. The effective rate k e fi 
of transition from / to III would then be given by The reflecting boundary was chosen 
to be at 9 = —100° where the relative free energy is ~ 10k bT. Using Eq. [5j we calculate 
the value of t(9i) for different pressures. We list the values of t(9i) in table [I] where D is 
measured in deg 2 / sec. We find that t(9i) increases upon increasing pressure within the error 
bars and increases sharply for P = 3000 atm. 
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Pressure 


t{ei)D- l l{f{deg 2 ) 


1 atm 


2.1625 


1000 atm 


9.8104 


2000 atm 


5.3807 


3000 atm 


1018 



TABLE I: Mean first passage times 



IV. HYDRATION SHELL AND BASE FLIPPING OF RNA 

As we have seen the sections above, that the base flipping kinetics changes as the pres- 
sure is increased - namely, the free energy barrier for the bulge base to flip out increases 
with pressure. Moreover, the transient barrier which presumably is due to the base triplet 
formation of the bulge base with the neighboring bases disappears at pressures P > 2 kbar. 
We note that 2 kbar is the pressure where most of the anomalies of liqid water disappears 
and also the pressure at which hydrophobic barriers for small molecules in water tend to 
vanish. Motivated by this we looked at the structure of the solvation shell (first hydration 
shell) of water around RNA for different pressures. We show a typical hydration shell around 
RNA in Fig. [61 The hydration shell is calculated by finding all the water molecules within 
a distance Rc of RNA molecule. We choose Rc = 0.223 nm as the first minimum in the 
radial distribution function of RNA and oxygen of water molecules (not shown here). We 
find that Rc is independent of the pressure. 

We first looked at whether the observed change in the pressure dependence of the kinetics 
is a result of ordering of water around RNA. To this effect, we calculated the OOO-bond 
angle and its distribution P(<f>) of a central water molecule and its nearest neighbors in 
the hydration shell. If the water orders then 0, would be very close to the tetrahedral angle 
109.47°. In Fig. [7J we show P{4>) for pressures P = 1, 1000, 2000, 3000 atm. For a comparison 
we also plot P(4>) for bulk TIP3P water at atmospheric pressures. We find that, the second 
peak corresponding to more ordered water of the distribution P(<p) shifts to smaller values 
of <j), suggesting that the water monotonically disorders upon increasing pressure. Moreover, 
we do not find any significant sharp changes in P(<p) which could be associated with sharp 
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FIG. 6: A typical hydration sheath around RNA. Water molecule is represented by a surface plots 
and the hydration shell is obtained as mentioned in section IV. 



change in the free energy barrier observed at P = 3000 atm. 

To quantify the orderin g /d isordering of water molecules around RNA, we use the tetra- 



hedral order parameter^ [24|-|26l|. Tetrahedral order parameter Q quantifies how close a 
given water molecule and its first shell neighbors form a structure close to a tetrahedron. In 
general, Qk of k th molecules is denned as 

Q" = 1 " I E X> os< ^ + V3) 2 (6) 

i>j j 

where, the indices i, and j run over all four neighboring molecules and 0^ the OOO-angle 
formed between the oxygens of central molecule k and neighbors i and j. 

Since, we only consider a thin hydration shell, the expression for ensemble average < Q > 
can be written as 

<Q>=l-^J (cos^ooo + l/3) 2 P(0ooo)#ooo (7) 

In table [Til we list the average tetrahedral order parameter < Q > of the hydration shell for 
different pressures. For a comparison, we also compute < Q > for bulk water. Table [TTJ lists 
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FIG. 7: Probability distribution function P ((f)) of the OOO-angle <j> of the hydration shell for 
various pressures. Note that a comparison with P(4>) for bulk water at P = 1 atm suggests, that 
hydration shell around RNA is more disordered and this disorder increases monotonically upon 
increasing pressure. 



Average Tetrahedral Order Parameter < Q > 


Pressure 


Bulk Water 


Hydration Shell 


1 atm 


0.467 


0.396 


lOOOatm 




0.375 


2000 atm 




0.361 


3000atm 




0.348 



TABLE II: Average tetrahedral order parameter Q of the first hydration shell of water around 
RNA at different pressures 



values of < Q > for different pressures. We find that, < Q > monotonically decreases upon 
and no sudden change in < Q > is seen at P = 3000 atm. 



increasing pressure 



Since, we did not see any sudden change in the ordering of hydration shell around RNA 
that might lead to the base flipping barrier observed at P = 3000 atm, we next studied the 
hydration shell of the bulge base and suprisingly, we find that at high pressures, the average 
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number of water molecules in the first hydration shell of the bulge base increases, from 
an average of about 1.0 to 1.40 (see Fig. EJ^b)). Moreover, we find that the distribution 
of water molecules around the bulge base shows significant probability of finding 3 — 4 
water molecules. To this end, we suggest that the increased barrier of base flipping and 
the disappearance of the transient barrier in the free energy barrier is due to the presence 
of increased water molecules in the solvation shell around the bulge base. The presence 
of more water molecules creates a competition between the base triplet formation and the 
hydrogen bond formation with the water molecules, which might have penetrated from the 
major groove, and hence stabilizing the stacking of the bulge base. 

U ' ~r 1 1 1 1 1 1 1 3 1.6i 1 1 1 1 1 



- P=latm 

- P=2000atm 




FIG. 8: (a) Probability distribution P{ny/) of water molecules in the first hydration shell around 
the bulge base "A" for different pressures. Note that the probability of 3 or 4 water molecules 
increases as the pressure is increased, (b) Ensemble averaged < nw > as a function of pressure. 



V. SUMMARY AND DISCUSSIONS 

In summary, in this paper we have investigated the effect of pressure on the kinetics of 
base flipping by calculating the free energy barrier associated with looping out torsional angle 
of an adenosine bulge base in a short double strand RNA. We find that, upon increasing the 
pressure, the propensity or likelihood of base flipping decreases. At pressure P > 2000 atm, 
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we see a sharp increase in the free energy barrier. Further, we calculate the time scale of 
flipping by mapping the problem of base flipping to a diffusion of reaction coordinate on 
an underlying free energy landscape from which we calculate the time scale of looping out 
of bulge base. We find that the time scale increases upon increasing pressure and changes 
dramatically at P > 2000 atm. We associate this behavior with the competing hydrogen 
bonds between the neighboring bases and increased solvation of the bulge base at high 
pressures. 
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